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Abstract 

We study the stochastic susceptible-infected-recovered (SIR) model with time-dependent forcing using analytic techniques which 
"^ "allow us to disentangle the interaction of stochasticity and external forcing. The model is formulated as a continuous time Markov 
^__l process, which is decomposed into a deterministic dynamics together with stochastic corrections, by using an expansion in inverse 
C^ system size. The forcing induces a limit cycle in the deterministic dynamics, and a complete analysis of the fluctuations about 
Cn| this time-dependent solution is given. This analysis is applied when the limit cycle is annual, and after a period-doubling when 
j^it is biennial. The comprehensive nature of our approach allows us to give a coherent picture of the dynamics which unifies past 
^ work, but which also provides a systematic method for predicting the periods of oscillations seen in whooping cough and measles 
■^ epidemics. 

■TJ" Keywords: non-linear dynamics, period doubling, measles 
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The availability of extensive time-series data for childhood 
diseases is often the reason given for the amount of atten- 
tion that this subject receives. However, the possibility that 
relatively simple models can capture the essence of the dis- 
ease dynamics also makes the topic an attractive one for mod- 
CN ellers. Mathematical epidemiologists are especially intrigued 
^ ] by the rich variet y of oscillatory dynam ics seen in this data 
^ llEarn et all I2OOOI: lOrenfell et all I2OO2I) . Wifliin the litera- 
^f-s Jure there is a broad consensus that there are two main ele- 
psj ments needed to model these oscillations: firstly st ochastic - 
ly-T ity, due t o the individual nat ure of the population (IBartlett , 
f^ [i960; Du rrett and Levinill994l) : and secondly, seasonal forcing, 
arising from the term-tim e aggregation of c h ildren i n schools, 
which is dete rministic (ILondon and Yorki Il973t ISchenzle . 
[1984ilAltizer e t al., 2006|). Independently these two factors are 
k> ^^^^ understood, but how they interact w hen both included in 
rS the same m odel is still an open question dKeeling et al.L uOOU 
iRohaniet al., 2002: ICoulson et al.ll2004l) . 

Measles is the canonical example of a disease which dis- 
plays recurrent epidemic behaviour. In larger cities regu- 
lar periodic oscillations, usually annual or biennial, are ob- 
served, whereas smal l er cities display more irregul ar dynamics 
dOrenfeU et al.l 120021 lUovd and Sattenspiell l2009h . The intro- 
duction of mass vaccination in the 1960s provides a 'natural ex- 
per iment' after which the dynam ics become much more irregu- 
lar (IGrenfell and Harwoodll 19971) . One of the early successes in 
the field was a simple deterministic model with seasonal forcing 
which could recreate the regular dynamics of measles dOietzi 
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1976tlSchwartz and Smithlll983b . Where simple models such 
as this fail, is in captu ring the more irregu lar dynamics seen 



in smaller populatio ns (IGrenfell et all 120021) . after vaccination 



(IRohani et all 1199^. and in other diseases such as whooping 



gh (" Nguven and RohaniLl2008h . These aspects can only be 
captured by fully stochastic models. 

Stochastic models without external forcing show large oscil- 
lations cau sed by the stochasticity exciting th e system's natural 
frequency ( lBartlettlll957l : lAlonso et al.ll2007l) . When forcing is 
included, it is less clear how the stochasticity interacts with the 
cyclic solutions that are produced. It could act passiv ely to kick 
the sy stem between different deterministic states (i Schwartz . 
19851) . as well as interacting with the non-hnearity to excite 



the transients. Power spectra dPriestlevl 1 1 982t [Anderson et al 



11984.) have proved a useful tool in investigating these factors. 
They can help distinguish various components in the time- 
series and classify them as essentially seasonal, stochastic or 
an interaction of the two (■BentorL.2006) . The most successful 
synthesis, by iBauch and EarnI d2003bh . showed that a simple 
mechanistic model can accurately predict the position of peaks 
in the power spectrum of a number of different disease time 
series. 

We approach this problem by starting with an individual 
based model, which is inherently stochastic. We can then both 
simulate it and derive the emergent population level dynam- 
ics. The novel aspect of this work is that we calculate the 
power spectrum analytically for explicitly time-dependent ex- 
ternal forcing, and compare the results with stochastic simula- 
tions. We do this by formulating the model as a m aster e quation 
which can then be studied using van Kampen's (1992^ expan- 
sion in the inverse system size. The macroscopic dynamics can 
then be viewed as a sum of a deterministic and a stochastic part. 
The value of the analytic approach is that we can more easily 
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deduce the mechanisms behind the dynamics and better under- 
stand the interplay between deterministic and stochastic forces. 
The theory we develop in this paper unifies much of the 
previous work on these models . It encompasses the influ- 
ential work of Earn et al.l (120001) in understan ding the transi- 
tions in measles epidemics, the later work of IBauch and EarnI 
(l2003bh relating to the transient fluctuations close to cycUc 
attractors for different diseases and the more recent work 



on stochastic amplification in ep idemic models dAlonso et al 
2007t Black and McKang^ 201 Of). The picture th at emerges is 



close to that proposed bv Bauch and EarnI (l2003bl) . but goes be 
yond it in two important respects. Firstly, we calculate the exact 
power spectrum for the forced model. Secondly, we show how 
the forcing changes the form of the fluctuations, and how in a 
stochastic model these are intimately related to the period dou- 
bling bifurcation, which is vital for explaining the dynamics of 
measles. 

The rest of this paper is as follows. In section |2] we intro- 
duce the seasonally-forced version of the stochastic susceptible- 
infected-recovered (SIR) model and the system-size expansion 
of the master equation. Section [3] provides a discussion of the 
results for the simple case when the deterministic dynamics are 
described by an annual limit cycle. In Section |4l we apply our 



method to elucidate the dynamics investigated by Earn et al 
(1200 0), which can account for the transitions in measles epi- 
demics. This is an interesting parameter regime, as the deter- 
ministic theory predicts a period doubling bifurcation. Finally, 
in section |5] a broad discussion of our results is given, describ- 
ing how this approach can account for the different dynamics of 
measles and whooping cough. There are two appendices giving 
technical details relating to the system-size expansion and Flo- 
quet theory. 

2. The seasonally forced SIR model 

We first summarise the individual-based stochastic SIR 
model. We emphasise only the aspects which are relevant to 
this paper; a m ore general discussion can be found in textbooks 
on the subject ( lAnderson and Mavi Il991t [Keeling and Rohani, 
12007) . The population is split into three classes: susceptibles, 
infected and recovered. Birth and death rates are set equal to 
yU and these events are linked, even in the stochastic model, so 
that the total population, A^, remains constant. Recovery hap- 
pens at a constant rate y, so that the average infectious time is 
1 /y; once recovered, an individual is immune for life. Seasonal 
forcing is included by assu ming that the tra nsmission rate Pit) 
follows a term-time pattern (ISchenzleLll984l) . 



y6(0=j6o(l+j6iterm(0), 



(1) 



where ySo is the basehne contact rate, [i\ the magnitude of forc- 
ing and term(f) is a periodic function which switches between 
1 during school terms and -1 during holidays. In this paper we 



use the England and Wales term dates set down bv lKeeling et al 
(1200 lb . The reproductive ratio is determined by Rq - ifi)ly. 
where {ji) is the effective (time-averaged) transmission rate: 



and ps is the proportion of time spent in school. For our choice 
of terms p^ - 0.75. We also include a small immigration term, 
77, to account for infectious imports. We use a commuter formu- 
lation, where susceptibles are in c ontact with a pool of infec- 
tives outside t he ma in population (lEngbert and Dreppen. 1 1994c 
Alonso et all |2007). Since A^ = 5 -t- / -H /?, we can use this 
constraint to eliminate the variable R from the rate equations. 

The model is then defined by the processes through which it 
evolves. If we write the state of the system as cr = {5,/}, we 
can specify the following transition rates, r(cr|cr'), between an 
initial state cr' and a final state cr: 



1 . Infection: 5 + / 



/?(') 



/ + / and 5 



r(5-l,/+l|5,/) = (A0-/ + '75) 



2. Recovery: I ^> R. 



T(SJ-l\SJ)^yI. 



3. Death of an infected individual: / — > 5 . 



T(S +l,I-l\SJ)^iuI. 



4. Death of a recovered individual: /? — » 5 . 



T(S + l,I\S,I)=fi(N-S -I). 



(3) 



(4) 



(5) 



(6) 



Since birth and death are coupled, the processes 3 and 4 also 
imply the birth of a susceptible individual. An important point 
is that changes in vaccinati on can be mapped onto the effective 
transmission rate (/?) using (lEarn et all l2000h 



(J3} ^ {fiXl - p), 



(7) 



where p is the proportion of individuals vaccinated at birth. We 
can also approximately map a change in birth rates onto (J3}, 
but this is not exact in this model because births and deaths 
are linked. In this paper we are primarily interested in the pa- 
rameter range of childhood diseases. These are characterised 
by yu «c 7, i.e. the average life expectancy of an individual 
is orders of magnitude lon ger than the mean infectious period 
(Anderson and May, 1991). 



2.1. Methods of analysis 

We use two methods to investigate the dynamics of this 
system . Firstly we simulate the system using Gillespie's 
(119761) algorithm with the appropriate time-dependent exten- 
sions (lAndersoni l2007h . This method generates exact realisa- 
tions from which statistical quantities, such as power spectra 
and moments, can be computed. The second method is ana- 
lytic, through the construction of a master equation. The master 
equation describes the evolution of the probability distribution 
of finding the system in state cr at time t. 



(J3}^MPs(l+/3i) + (l-Ps)(l-/3i)l 



(2) 



dPio-; t) 
dt 



= ^ T{cr\cr')P(a'- f) - ^ T {a' \a)P(cr- 1). (8) 



a'i^cr 



This c annot be solved exactly so we instead use van Kampen's 
(Il992h expansion in the inverse system size to derive approxi- 
mate analytic solutions. This involves making the substitutions, 



as the system alterna t ely sw itching between two spiral fixed- 



S ^N<p + N"^x, 
I^Nip + N^'^y, 



(9) 



and expanding the master equation in powers of A^ '''^. This 
technique and similar ones have been documented at length 
in the li terature, but almost exclusively for time-independent 
models (lAlonso et all l2007h . The novel aspect o f this paper is 



that w e analyse the full time-dependent system jBoland et al. 



20091) . whereas in a previo us paper we used the appro ximation 
which replaced y6(0 by (/S) dBlack and McKanellloiol) . 

The details of the system-size expansion for this model are 
given in Appendix A. At leading order we find a pair of deter- 
ministic equations, describing the mean behaviour, which scale 
with the system size A^, 






(10) 



These are the same as the equations that are found from a purely 
phenomenological treatment of the SIR system. At next-to- 
leading order we obtain a pair of Langevin equations for the 
stochastic corrections to the deterministic equations (fTOl i. 



X = K(t)ii{t) + f(0, 



(11) 



where x = {x, y), and f (f) are Gaussian white-noise terms with 
correlation function {i{t)t{t'Y) - G{t)6{t - t'). The matrices 
K(t) and G(t) are determined from carrying out the expansion 
and are given by 



K(t) 



-pip- T]-fl 

P^ + ri 



■y-H 



(12) 



and 



Gil =yS<;*iA + # + Ml -0). 
G22=j8# + # + (r + /")'A, (13) 

Gl2 = G21 = -y8# - # - A"A> 

where a bar indicates that the solutions are evaluated on the 
limit cycle. These are essentially the s ame as are found in the 
non-forced model dAlonso et al.L l2007b . except now /?, and 1^ 
are all functions of time. 



3. Stochastic amplification about a limit cycle 

The mean behaviour is found by integrating the deterministic 
equations (fTOll. WhenySi = 0, solutions show damped o scilla- 
tions tending to a fixed point (lAnderson and Mavi Il99lh . For 



non-zero /3i , t his model can disp l ay a rich set of dynamic s in- 
cluding chaos dOlsen et al.[ll988l:lRand and Wilsonlll99lh . but 
for realistic parameter values the most common long-time so- 
lution is a limit cycle with a period that is an integer m ulti- 
ple, n, of a year dOietzl Il976t ISchwartz and Smithl Il983h . As 
the forcing is a step function in time, we can visualise this 



points (iKeeling et al.L 1200 lb resulting in a piecewise contin- 
uous limit cycle, illustrated in figure [1] Any other periodic 
forcing function, for instance a sinusoidally varying one, could 
be used without more difficulty, and would typically lead to a 
limit cycle which is smooth. As /3i is increased, the limit cy- 
cle grows (although typically not linearly with /3i) and at crit- 
ical values bifurcations are induced to longer period solu tions 



dAron and Schwartzlll984tlKuznetsov and Piccardill 19941) . 
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Figure 1 : Phase portrait illustrating a deterministic solution of the forced SIR 
model. The term-time forcing creates a limit cycle (red curve) as the system 
alternately spirals between the two fixed points defined by /J// =/?o(l +/?i) and 
jSl = /So(l ~/?i )■ The light blue solutions show the behaviour if the forcing was 
switched oft', to illustrate the two spiral attractors. The red dot shows the fixed 
point calculated using the approximation where jS(;) is replaced by (/}}. 

In this section we present results where there is only an an- 
nual limit-cycle (n = 1). The case where we also have a pe- 
riod doubhng is examined in section H) The stability of these 
limit-cycle solutions can be investigated with the use of Floquet 
theory (jKuznetsov, 2004; Boland et al., 2009). This quantifies 
how perturbations to the trajectory of the limit cycle behave 
and is analogous to linear stability analysis about a fixed point 
(Grimshaw, 1990). 

Floquet theory states that for any periodic solution of Eq. ( flOb 
there exists a matrix B which satisfies the relation. 



X(t + T„) = X{t)B, 



(14) 



where X{i) is the fundamental matrix (IGrimshawl Il990h and 
T„ is the period of the limit cycle. The eigenvalues of B are 
called the Floquet multipliers, p,; a related set of quantities are 
the Floquet exponents /i, = ln(p,)/r„ (in this paper, since we 
will be discussing frequencies rather than angular frequencies, 
these exponents will be divided by a factor of 2;7r). Another way 
to think of this is as linear stability analysis o f the fixed points 
of the n-cycle Poincare map of the system (B auch and Earn , 
|2003b; Kuznetsov, 2004) . A limit-cycle solution will be stable 
if \pi\ < 1. When the multipliers are complex, perturbations 
to the trajectories return to the limit-cycle in a damped oscilla- 
tory m anner, analogous to a stable spiral fixed point (GrimshawL 
1990h . Similar ideas have been used to investigate the transients 



in forced epidemic systems in the past, but only in a determinis- 
tic setting (Bauch and Earn, 2003b ; He and Earn, 2007) . Here 
we wiU explore how the nature of the fluctuations can be quan- 
tified using Floquet theory. 

Figure |2] shows a simulation of the full stochastic system 
together with the deterministic limit-cycle solution. We can 
see that even for large populations the stochastic corrections 
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Figure 2: Phase portrait of the stochastic SIR system. A time-series of 100 
years duration is shown in light blue. The first two years are highlighted in 
dark blue, with the dot showing the start point. The macroscopic limit-cycle 
(red) is also superimposed Parameters are those relevant for whooping cough 
jNguven and Rohanil20og) : i?o = 17, y = 1/22, /3i = 0.25, ^ = 5.5 X lO"', 
n = lO"*^ and Af = 2x10*. 
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Figure 3: Power spectra for the number of infectives from simulation (Ught 
blue solid curve) and analytic calculation (black dashed curve). From the sim- 
ulations, we observe a sharp peak at 1 year from the deterministic annual-limit 
cycle. The other peaks, marked by the red lines, are from stochastic amplifica- 
tion, with the peak frequencies given by m ± Im(/i), where lm{A) = 0.36. The 
dominant stochastic period is therefore 1/0.36 = 2.7 years. Parameters ai'e as 
in figure 12] 



to the deterministic solution are important. The noise due to 
demographic sto chasticity (noise at the in dividual level due 
to chance events; iNisbet and Gurnevll 19821) excites the natural 
oscillatory modes about the limit cycle, creating a resonance 
and giving rise to large scale cohere nt oscillations — an ef- 
fect known as stochastic a mplification (IMcKane and Newman . 
2005t lAlonso et all l2007h . As described in Appendix B, by 
solving Eq. dTlT ). and invoking aspects of Floquet theory, we 
can express the auto-correlation function. 



C(t) = - 

^ n Jo 



(x(f + t)x' (0) dt, X = {x, y]. 



(15) 



as an integral without further approximation ( IBoland et al. , 
20091). Taking the Fourier transform of this expression then 
gives an exact expression for power spectrum of these stochas- 
tic oscillations. 

Figure |3] shows simulated and analytic power spectra for the 
system shown in figure |2] We observe a sharp peak at 1 year 
due to the deterministic annual limit-cycle and a number of 
broader peaks due to the stochastic amplification of the tran- 
sients. We would expect on general grounds that the stochastic 
peaks would be observed at frequencies. 



m/T„ + lm(A), 



(16) 



where m is an integer and A is the Floquet exponent JWisenfeld , 
1985b : Rolan d et al.ll2009h . and this is indeed what is seen. For 
the annual limit-cycle the dominant peak is at + lm(A), with 
the other peaks being much smaller. Near to bifurcations these 
minor peaks become important and are treated in more detail in 
the following section. 

The area under the peaks in the power spectrum is propor- 
tional to the root-mean-square amplitude of the oscillations. 
Away from any deterministic bifurcation points the amplitude is 
proportional to Re(/i), as in the unforced model. Thus the spec- 
trum is close in form to that predicted from the unforced model 



by substituting (B) for the time-independent transm ission rate 
dBauch and Earn[l2003at black and McKanell2010h . 

There is good agreement between analytical calculations and 
simulations. Although calculations give the power spectrum 
as an integral, it must be evaluated numerically because the 
deterministic equations ( fTOb cannot be solved in closed form; 
this is all carried out using th e symbolic package Mathemat- 
ica (I Wolfram ResearchL 120081) . T liis analysis about an an nual 
limit cycle corresponds to that of iBauch and EarnI (l2003bl) ex- 
cept that we can derive the full power spectrum. They term the 
'resonant peak' what we describe as the deterministic or an- 
nual peak, and the 'non-resonant peak' what we describe as the 
stochastic peaks. Their terminology is somewhat misleading, 
as the stochastic peak is generated by a resonance phenomena 
whereas the macroscopic peak is not. 



4. Period doubling and measles transitions 

We can use our analytic methods to help understand the 
dynamics and large-scale temporal tr ansitions in measl es epi- 
demic patterns, first investigated by ' Earn et alj (12000). The 
main force in driving these transitions is changes in the sus- 
ceptible recruitment (a mixture of changes in birth rates and 
vaccination), which can be mapped onto Rq. Thus a knowledge 
of the model dynamics as a function of Rq can be used to ex- 
pl ain the changes in epidemic patterns. Although the analysis 
of lEarn et alj ( l2000b is in good qualitative agreement with time- 
series data, there are a number of outstanding questions with 
regard to the interpretation of the mechanisms for the dynam- 
ics. We first provide a brief review of the original analysis and 
then go on to show how the stochastic dynamics of this model 
can be understood within the framework we have laid out in the 
previous section. 




Figure 4: Bifurcation diagram sliowing tlie SIR dynamics as a function of i?o. 
Fixed parameters: ;Si = 0.29, y = 1/13, ^ = 5.5 X 10"' and ;; = 0. The 
different period limit cycles are shown in different colours, which are produced 
by different initial conditions. 



As ri is decreased further still, more of t he structure is found 
(iBolker and GrenfellLll995l:lNasellLl2002h . 

Immigration is an important aspect in the simulation because 
without it the disease would fade out as the minimum number of 



infections can go far below a single individual (IBartletti 11957 ; 



BolkerandGrenfell. 1995: Conlanetal.. 2009). In a determin 



istic analysis this term is easily omitted because the v ariables 
are c ontinuous and therefore fadeout cannot happen dNasell , 
19991) . This raises the question: do these longer period solu- 



tions have an effect on the stochastic dynamics? If not, how can 
we describe the nature of the stochastic dynamics? We can use 
our analytic method to help clarify these questions. The power 
spectrum is especially useful as it can show up anomalous peaks 
from simulations. 



4.1. Review of original analysis 

It is acknowledged that stochasticity plays a role in the dy- 
namics of measles, which can only be captured through sim- 
ulation of the individual-based mod el. Fundarnentally though, 
the analysis of these mechanisms bv lEarn et al.l (120001) is deter- 
ministic. Figure m shows the bifurcation diagram derived from 
the SIR equations ( fTOt . as a function of /?o, with parameters 
corresponding to measles and no immigration (77 - 0). This 
shows the incidence sampled annually on the 1st of January 
each year, thus stable limit cycles are shown by different num- 
bers of (colour coded) curves. The single curve, beginning at 
small /?(), shows an annual cycle which bifurcates at /?(> = 15.5 
into two curves giving a biennial cycle. For values of 7?o ly- 
ing between about 5 and 15, there are several sets of n curves 
representing n-year cycles. 

For large Rq (e.g. Rq x; 30) only an annual limit cycle exists. 
As Ro is reduced a biennial limit-cycle is found; before vacci- 
nation was introduced in England and Wales, most cities would 
be in this region. Higher birth-rates might move the system 
back into the region with only an annual attractor, whereas vac- 
cination would act to reduce Rq, moving it into the region with 
multiple co-existing longer period attractors. The interpretation 
put forward by Earn et al, is that stochasticity will then cause 
the system to jum p between these different deterministic states 
(ISchwartall985h . giving rise to irregul ar patterns. Thus, in this 
description, noise plays a passive role dCoulson et al.1 120041) . 

Although peaks were seen in power spectra from simula- 
tions, which appear to confirm this view, there are a num- 
ber of problems with this interpretation. The crucial aspect 
that is neglected is that there are no infectious imports in- 
cluded in the deterministic analysis (although presumably they 
are included in simulations). When this factor is introduced 
(77 ^ 0) then most of the additional structure disappears, see 
figure |5^; we are left with an a n nual hmit cycle and a period 



4.2. Analytic predictions 



1996; 



doubling (Enabert and Dreppen 119941 iFerguson et al 
Alons o et al., 2007). 

When T] - 10"^, there is only a small region in the range 24 < 
Ro < 25 where there are coexisting annual and biennial limit- 
cycles. As the immigration parameter is reduced some of the 
additional structure reappears; for example at 77 = 10"^ some of 
the period 3 attractors can be found in the range 9 < Rq < 11. 




Figure 5: (a) Bifurcation diagram for the SIR model with fS[ = 0.29 and 
;; = 10"*. (b) The Floquet multipliers, which are in general a complex con- 
jugate pair, thus we plot the real (dark blue) and imaginary (light green) parts 
separately, (c) Imaginary parts of the Floquet exponents. Note that in the region 
where there are the coexisting limit cycles (24 < Rq < 25), only the multipli- 
ers/exponents for the biennial cycle are shown for clarity. 

Figure |5] shows the bifurcation diagram for the model pre- 
sented in the previous section, but with rj - 10 '', along with 
the Floquet multipliers and exponents. These parameter values 
will be used for the rest of this section. Figure |6] shows the 
Floquet multipliers on a larger scale near the period doubling 
bifurcation point and figure |7] shows the analytical and numer- 
ical power spectra for various values of 7?o with N = 5 x 10^. 
Away from any bifurcation points there is good agreement be- 
tween the analytic and the simulated spectra. 

As we approach the period-doubling bifurcation point from 
belo w, the stochastic o scillations follow a virtual-Hopf pat- 
tern JWisenfeldl ll985alb V This is where the oscillations first 
show the precursor characteristics of a Hopf bifurcation be- 
fore changing into the precursor characteristics of a period- 
doubling. This is clearly seen in the power spectra shown in 
figure|2] In the Hopf-like regime (/?() < 14.94), the Floquet mul- 
tipliers are a complex conjugate pair, giving rise to two peaks 
in the spectrum: a major one at frequency Im(/l) and a minor 




(a) 



Figure 6: Floquet multipliers near to the period doubling bifurcation point, 
showing the virtual-Hopf pattern. For Rg < 14.94 the multipliers are a complex 
conjugate pair, with a negative real part (dark blue line); this is the Hopf-like 
region. The actual period-doubling bifurcation occurs at RJJ'' = 15.34, where 
one of the multipliers becomes equal to - 1 . 



one at l-Im(/l), as in section [3] Therefore in figure |7^ the two 
peaks are most widely separated for Rq = 4. 

As we increase Rq, lm(A) also increases, and the major and 
minor peaks move closer together, converging at 0.5 y"' when 
the multipliers become real and negative; this marks the on- 
set of the period doubling regime, see figure |6] In this regime 
(14.94 < Rq < 15.34), as the multipliers are negative, so their 
phase is +n and so the imaginary part of the Floquet expo- 
nents is +nl2jiT\ - ±0.5. Therefore the peak stays fixed at 
0.5 years"' as we increase Ro further within this range, but the 
amplitude increases quickly. At R^'^ - 15.34 one of the mul- 
tipliers reaches -1 and we see a deterministic period doubling 
JKuznetso v. 2004), and the size of the fluctuations grows to or- 
der A^. Figure[8]shows how in this way the oscillations smoothly 
turn into the macroscopic biennial limit cycle. The same pattern 
is seen if we hold Rq fixed and increase /3i to induce a period 
doubling. 

When the system is in the biennial regime we can still cal- 
culate the fluctuations about the limit cycle and get a good 
correspondence with analytic predictions (figure [TJ)). The po- 
sitions of the peaks are now at m/2 + lm{A) and the spec- 
trum changes little within this parameter range. The peaks at 
m/2 + lm{A) are barely visible, as compared to the prominent 
peaks at m/2 - lm(A). In the annual regime after the doubling 
(Ro > 25), the analytic results are again very accurate, with 
stochastic peaks at frequencies m + lm(A) (figure |7};). Here as 
well, the set of peaks atm + Im{A) are much smaller. Note that 
in both of these regions the time-series will be dominated by 
the deterministic signal as the stochastic oscillations are much 
smaller than in the pre-bifurcation region {Rq < 15). 

4.3. Near the bifurcation point 



For values of Rq near the bifurcation point, the deviations 
between the analytic and simulated spectra become larger (see 
for example figure |7h; ^o = 14). This is expected: the anal- 
ysis developed here is essentially linear and thus predicts an 
unbounded incr ease in the flu ctuations as we approach the bi- 
furcation point dGreenman and Benton. 2005) . As the fluctua- 
tions become larger the Unear approximation breaks down and 
non-linear effects become important and act to bound the fluc- 
tuations. Going to larger system sizes can result in better agree- 
ment between analytic results and simulation, but this will al- 
ways break down at some point. 
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Figure 7: Analytic (black dashed curves) and numerical (coloured) power spec- 
tra for a range of Rq with W = 5 X 10^. In most cases the analytic and numerical 
spectra are virtually indistinguishable, apart from Rq = 14. (a) Spectra before 
the bifurcation, Rq = 4,6,8, 10, 12, 14. (b) Typical biennial regime, Rq = 20. 
Note that the stochastic peaks have been made clearer by subtracting the deter- 
ministic dynamics before calculating the power spectrum. The spectrum would 
otherwise be dominated by the peak at 0.5 y"' . (c) The major and minor peaks 
in the large Rq annual regime: Rq = 26, 30, with the larger peaks corresponding 
to Rq = 26 for both the major and minor peaks. 



Although the analytic approximation breaks down near the 
bifurcation point, the structure we have uncovered is still vis- 
ible. Figure [8] shows stochastic power spectra from simu- 
lations for 14 < Rq < 18, as we move though the bifur- 
cation point. The virtual-Hopf pattern is still clear, as pre- 
dicted by the analysis, but the fluctuations re main bounded , 
growing to the same order a s the system size (van Kampenl 
1992: Kravtsov and Surovvatkina. 2003). Within this region the 
macroscopic dynamics cannot be split into a deterministic and 
stochastic part and it is not in general possible to reconstruct the 



deterministic part by averaging over many realisations. Thus, 
dete rmining exactly wh ere the bifurcation takes place is diffi- 
cult (IWisenfeldlll985bl) . At Rq - 16 the deterministic biennial 
peak should be observed, but is not clearly visible until Rq - 18. 
It is possible that the bifurcation point is shifted in the stochas- 
tic system, but more analysis is required to determine that this 
is so. 
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Figure 8: Simulation results showing how the power spectrum of the stochastic 
oscillations changes as the period doubling bifurcation point is crossed. The 
peaks for Rq = 17 and 18 have been cropped for clarity. 



4.4. Smaller populations 

The results presented in the previous sections were for A^ = 
5x10^, which roughly corresponds to the largest populations 
we would be interested in modelling. Simulations of smaller 
populations tend to show regular deviations from the analytic 
calculations and results are sensitive to A^, 77 and ySi . The forc- 
ing pushes the system close to the fade-out boundary (/ - 0), 
where fluctuations are non-Gaussian, and so large deviations 
from the theory are expected. Figure |9] shows the stochastic 
power spectra from simulations, within the range 4 < Rq < 30 
and with A^ = 5 x 10^ 10^ and 5 x 10\ 

For smaller values of Rq we still clearly observe the virtual- 
Hopf pattern, but a visual inspection of the time-series shows 
much more irregular dynamics. This is due to the increased 
stochasticity in the smaller systems, but also the closeness of the 
fade-out boundary, where extinction and re-colonisation events 
start to have an impact on the dynamics (Gr iffiths, 1973) . This 
has an efifect on the power spectra in two ways: firstly as a 
broadening of the power spectra, showing a greater range or 
amplified frequencies and thus a more irregular dynamics. Sec- 
ondly the endogenous period is syste matically shifted hig her, as 
in unforced versions of this model dSimoes et al.l uOOSh . This 



reflects the fact that the period of oscillations also depends o n 
the re-introduction of the disease after fade-out (lBartlettlll957h . 

The most important effect is on the fluctuations in the bien- 
nial regime after the period doubling. For A^ = 5 x 10^ the 
peaks are sharp, indicating a deterministic limit cycle, and the 
stochastic oscillations are much smaller (figure Uh), hence the 
good predictability of these larger systems. For the two smaller 
populations this is not the case. We do not observe the de- 
terministic biennial limit cycle, but instead see an enhanced 
stochastic peak and a broadening of the spectrum. The range 
of this enhanced region is also reduced. 

Although there are large deviations, having an analytical de- 
scription still helps us interpret the dynamics at smaller A^. Tak- 
ing the average of Eq. (|9]l we obtain 



{I)^Nm+NHy). 



(17) 



In the linear noise approximation, which we have used in this 
paper, the fluctuations are Gaussian and therefore (y) - 0. At 
some point this will break down and we must include the next 
order corrections, which will be of the order A^'^^, to the macro- 
scopic equatio ns. It will no longer be true that th e mean is equal 
to the average dvan Kampenlll992tlGrimail2009h . This effect of 
the fluctuations on the deterministic dynamics could be enough 
to retard the onset of the biennial limit cycle and is the subject 
of further research. 

4.5. Switching between attractors 

As seen from the bifurcation diagram in figure |5] where 
77 = 10"^, the only region where deterministically there are pre- 
dicted to be two coexisting states is when 24 < Rq < 25. This 
can be detected in simulations and the period of switching de- 
pends strongly on the system size. If the system is large it will 
tend to stay in the state it started in, because the fluctuations are 
not large enough compared to the mean to kick the system into 
the other state. Decreasing the system size makes this possible, 
and we see periods of annual dynamics followed by biennial 
and back to annual, where the period of switching depends on 
the system size. 

There is another intriguing region where we see signs of this 
type of behaviour. For A^ = 10^ and Rq = \0 (figure |9]5), we 
observe an enhanced stochastic peak in the spectrum with a pe- 
riod of 3 years. Visual inspection of the time-series shows re- 
gions of irregular annual oscillations interspersed with very reg- 
ular triennial oscillations. Note that this is not observed in the 
larger or smaller systems and the power spectrum is shifted by 
the proximity to the fade-out boundary from its infinite system 
size limit. Very similar behaviour is observ ed for measles data 
from Baltimore betwee n 1928 and 1935 (ILondon and York , 



'1973; Earn et al. 



(iBauch and EarnLl2003bi) . 



J2OOO), which has similar parameter values 



5. Discussion 

We have used an analytic approach, as well as simulations, 
to quantify the effect of stochastic amplification in the forced 
SIR model. The time dependence has been treated explicitly. 
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Figure 9: Power spectrum through the bifurcation point for different size populations, (a) N = 5x 10*, (b) N = 10*, (c) N = 5x 10^. Some of the peaks are cropped 
for clarity. Notice the anomalously enhanced peak for N = 10*, Rq = 10, see section l431 for discussion of this. 



instead of by using an app roximation as in a previous paper 
(IB lack and McKaneL l2010l) . Because this model has a finite 



namics show a strong depe ndence on population size jBartlett , 



1957; Grenfell et al. 



2002h . Our analysis also offers some in- 



population, and therefore is inherently stochastic, it can only 
be studied 'exactly' by simulation. The system-size expansion, 
which we use to derive approximate analytic solutions for this 
model suggests that we should view the population-level dy- 
namics as being composed of a deterministic part and a stochas- 
tic part, where the spectrum of the stochastic fluctuations is inti- 
mately related to the stability of the deterministic level dynam- 
ics. Power spectra of these models have been known for some 
time, but it has not always been clear what the mechanisms that 
generate the peaks are. This is the main advantage of being 
able to calculate the power spectrum of the stochastic fluctua- 
tions analytically; by comparison with the simulations we can 
gain insight into the mechanisms at work. 

Our analysis suggests a simple explanation for the dififer- 
ences seen in the epidemic patterns of measles and whooping 
cou gh in England an d Wales both before and after vaccination 
(iRohani et al.L 1 1999b . and which are representative of the two 
main parameter regimes for childhood diseases. The generic 
situation occurs when we are far away from the bifurcation 
point. Here we observe a deterministic annual limit cycle with 
stochastic oscillations, as in Figures|2]and|3] In general the form 
of the spectrum is close to that predicted by the unforced model. 
As already shown in a previous paper, this situation can account 
for the dynamics of whoop ing cough pre- and post-vaccination 
(IBlack and McKanel I2OIOI) . Pre-vaccination the stochastic os- 
cillations are centred on 2-3 years. Vaccination acts to shift 
the endogenous frequency lower and increases the amplitude of 
these fluctuations giving large four yearly outbreaks. 

Measles epidemics show a contrasting behaviour and rep- 
resent the second important parameter regime, where the de- 
terministic dynamics are near to a bifurcation point. Pre- 
vaccination, large cities such as London are in the regime with 
a deterministic biennial limit-cycle. Vaccination acts to lower 
Ro and shift the system into the regime where there is an an- 
nual limit cycle with large stochastic oscillations. As vaccina- 
tion coverage is increased, the endogenous perio d of these os- 
cillations is also increased dGrenfell et al.l 1200 ih . Measles dy- 



sight into this: in large populations the stochastic oscillations 
are very small compared with the deterministic biennial limit- 
cycle. This accounts for the regularity and explai ns why purel y 
deterministic models capture this aspect so well ( lDietzill976h . 
For smaller populations the deterministic biennial limit-cycle 
is not observed, just enhanced stochastic oscillations, thus ac- 
counting for the more irregular dynamics seen in these smaller 
populations. 

Finite size efi'ects and immigration / imports are closely re- 
lated in a stochastic individual-based model because the popu- 
lation is finite. Immigration reflects the basic fact that no pop- 
ulation is isolated and there must be reintroduction of the dis- 
ease if it fades out. One advantage of the approach which starts 
from an individual based model and derives the population level 
model, is that the immigration terms from the stochastic model 
are automatically included in the deterministic equations on the 
macro-scale. As we have shown, this is vital, as the longer pe- 
riod solutions are no longer found when they are included, and 
thus removes some speculation as to the influence of multiple 
co-existing attractors. These terms are easily omitted in a de- 
terministi c analysis be cause the system size is an innocent pa- 
rameter (lNaselllll999h . It is interesting to note the similarities 
between immigr ation and age-structure in these forced models 
( ISchenzleill984l) . Adding age-structure creates a constant pool 
of i nfectives in the infant class which acts to damp the dynam- 
ics (iBolker and Grenfell, 1993), exactly analogous to how we 
model immigration. 

Owing to the importance of immigration in these epidemic 
systems one relevant outstanding question is: what form should 
the immigration parameter take? As measles dynamics can 
be highly synchronised it could be argue d that the immigra- 
tion parameter should reflect this (Bolker and Grenfell. 1995t: 
Xia et all 120041: iLlovd and SattenspielL l2009l) . On flie other 



hand for larger cities, this parameter can be viewed as an 
aggregate of many infectious encounters from varied sources 
and could be approximated as a constant. Pre vious work has 
hinted at the sort of spatial effects that can arise (lKeelingll2000l: 



Keeling and RohaniLl2002tlHagenaars et all 120041) . but investi- 
gation of explicitly spatial stochastic models should be a high 
priority. 

The bifurcation diagrams for SEIR and SIR models are 
very similar, which justifies our use of the SIR model in 
this paper The extension to uncoupled births and deaths 
would be straightforward, but would offer no further insight 
dKeeling and Rohanil 120070 . There are technical difficulties in 
extending the method to the SEIR model because of the differ- 
ence in time scales between the collapse onto the centre mani- 
fold (ISchwartz and Smitlull983h and the period of forcing; this 
creates difficulties in computing the Floquet multipliers. These 
could in principle be overco me either by calculating the multi 



1991; Lust, 



pliers by a different method dFairgrieve and Jepson , 

20011) , or by carrying out a centr e-manifold reduction b efore 



doing the van-Kampen expansion (IForgoston et al.Ll2009l) . The 
breakdown of the linear theory near the bifurcation point can be 
remedied by including next-to-l eading order terms f rom the ex- 
pansion of the master equation dvan Kampenlll992h , but would 
result in a much more complex calculation. The calculations 
and discussions we have given here once again highlight the 
important role that simple models play in the understanding of 
complex systems. It also makes a novel contribution to the 
wider debate on the relative importance of stochastic and de- 



terministic forces in ecology tCoulson et am2004.) 
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Appendix A. Expansion of the master equation 

In this appendix we review the van Kampen system-size ex- 
pansion, which approximates the master equation ^ by the 
deterministic equations ( fTOl l plus stochastic fluctuations given 
by Eg. ([TTI) when A^ is large. The method has been described 
bviAlons o et al.l (l2007h for the SIR model without forcing (i.e. 
where /3 is independent of time), where further details are given. 
We begin by introducing step operators which allow us to 
express the master equation in a more compact form and also 
allow us to carry out the expansion in a more straightforward 
way. These are defined by: 



Ef\fiS,I,t) = f(S,I+l,t), 
Eff(S,I,t)^f(S±lJ,t). 

The master equation ^ can then be written in full as, 

+ (E-s'-l)[fi(N-S-P,] 



(A.l) 



+ (EiE-/ - 1) [///] + {Ej-l) [yl] \P(S,I, t), 

(A.2) 

which agrees with lAlonso etal.l (120071) . up to typographical er- 
rors in that paper. The essential step in the expansion is the 



ansatz (|9]); we anticipate that the approximate probability dis- 
tribution has a mean which scales as A^ and width which scales 



as N^^^ (I van Kampenll 19921) . To expand Eq. ( IA.2I ) in a power 
series in N^^^^, we write the step operators in terms of the fluc- 
tuation variables x and y: 



Ef = 1 +N'-2— + ^N' 



Ef' 



1 +A^" 



dx 
.d_ 

' dy 



kN- 



dy^ 



(A.3) 



Substituting these and the ansatz (|9]) into Eq. ( IA.2I ) we identify a 
hierarchy of equations multiplied by different powers of N^^^^. 
At leading order we find the deterministic equations (fTOl i for the 
macroscopic variables, (p(t) and i//(t). At next-to-leading order 
we obtain a linear Fokker-Planck equation for the fluctuations 
variables x s {x, y], of the form. 



f-Z-'/'> 



i,j 



d[xjU] 
dxi 



^Z^'vw 



'J 



dm 

dxidxj ' 



(A.4) 



where the matrices K{t) and G(t) depend on time through y6(f), 
(p{t) and ifr{t). The Fokker-Planck equation ( IA.4I ) is equiv- 
alent to the Langevin equations ( fTTT l given in the main text 
(van Kampen. 1992; Gardiner. 2003) . Since we will not be in- 
terested in transient behaviour, only fluctuations about the limit 
cycle, the solutions to ( fTOb will be those of the limit cycle, 
which we denote by (f> and ifr. The explicit forms for K(t) and 
G(t) are then given by Eqs. ( fT2] ) and ( fT3l l respectively. Since 
/3{t), (f){t) and i^(f) are periodic, so are K{t) and G{t). 

Appendix B. Floquet analysis 



The analysis of the Langevin equation ( fTTI ) is not difficult when 
the deterministic system approaches a fixe d point at large times , 
and for the SIR model this is described bv lAlonso et al.l (l2007h . 
When the attractor of the dynamics is a limit cycle, the analy- 
sis is more complicated, but can still be carried out to give the 
power spectra as integrals over kno wn functions. Here we out- 
line this analysis, following closelv lBoland et alJ (l2009h . which 
may be consulted for further details. 

We begin by considering linear perturbations about the limit 
cycle, that is, Eq. dTTT l. but without the noise which originates 
from the discreteness of the individuals. The equation describ- 
ing these small perturbations x = {x, y] is 



X = K(t)x, 



(B.l) 



where K(t) is given by Eq. (fT2l) . A fundamental matrix, X(t), 
is constructed from the linearly independent solutions of the 
homogeneous equation dB.lb . thus it satisfies the relation. 



X(t) = K(t)X(t). 



(B.2) 



The matrix X{t) is not unique and will depend on the initial con- 
ditions. Floquet theory states that if K(t + T„) - Kit), then there 
exists a canonical fundamental matrix which can be expressed 



in the form Xo(f) = P(t)Y(t) (lGrimshawlll990h . It has the prop- 
erty that the Floquet matrix, defined by Bq = XQ^(Q)Xo(Tn), 
is diagonal. These diagonal elements are called Floquet mul- 
tipliers. Pi and play a central part in the theory. The matrix 
P(t) carries the periodicity of the limit cycle, while Y(t) = 
diag[e'^'',e''-'], where A/ are the Floquet exponents. Using the 
periodicity of P{t) it can be seen that B = Y'^(Q)Y{T„) and so 
the Floquet multipliers are related to the Floquet exponents by 
Pi = e^'^",i = 1,2. Using the canonical form, we can derive 
an expression for the power spectrum in terms of the matrices, 
P(t) and Y(t). 

In practice, one obtains a fundamental matrix, X(t) by numer- 
ically integrating ( IB.2l i with initial condition X(Q) = I. With 
this choice of initial condition, B - X{T„). The multipliers 
and exponents are then calculated from the the eigenvalues of 
B, which allows the construction of the matrix Y{t), since the 
eigenvalues are independent of the choice of fundamental ma- 
trix. One can then calculate the canonical form, Xo(t) = X(t)S , 
where the columns of S are the eigenvectors of B. Finally P(t) 
is found from P(t) = X(0)Y-\t) = Xo(t)Y(-t). 

Having described the basic idea behind Floquet theory, we 
can now return to the Langevin equation ( fTTT l. which is an inho- 
mogeneous linear equation with periodic coefficients. We can 
use Floquet theory to construct a solution to this by adding a 
particular solution to the gene ral solution of the corresponding 
homogeneous equation ( IB. Il l dGrimshawl 1 1 9901) . This gives. 



xit) ^ X(t)X-\to)xQ + X(t) X-\s)fis)ds, 



J to 



(B.3) 



with initial condition x(fo) - xq. We are interested in the steady 
state solutions, when transients have damped down, thus we 
can ignore the first part of Eq. ( IB.3l l and set the initial time 
to the infinite past, fo — ^ -oo. Taking the case where X(t) is 
X()(t) - P(t)Y(t), one finds using the properties of the diagonal 
matrix Y(t), that 



kJ-O 



x(t) = Pit) Yit - s)Pisy'f(s)ds. 



(B.4) 



The correlation matrix is defined as C(t + T,t)- {x{t + T)x^(f)), 
which using Eq. ( IB.4b may be written as 

Xt+T r-t 
I Y{t + T-s)P(sr^G{s) 
OO <J —oo 

X 6(s - s')iP-\s')fY(t - sY ds' ds P(t)'^ , 

(B.5) 

where {t(s)f^(s')) = G(s)6(s - s'). Integrating over the delta 
function, the result will depend on the sign of r. If we take 
T > then the integration region is -oo < s < t, giving 



•f 

•J —OO 



C{t + T,t)^ P{t + t) \ Y{t + T- s)T(s)Y(t - sf ds P(tf , 



where we have defined 



r(.) = p{s)-'G{s){p-Hs)f , 



(B.6) 

(B.7) 



which will have the periodicity of the limit cycle. Next we make 
a change of variables, s ^ t - s' , which gives 



C{t + T,t)^P{t + T) Y(t + s')T{t - s')Y{s'f ds' P{tf . 
Jo 

(B.8) 

The form of Y means we may write Y{t + s') - Y{t)Y{s'), and 

so the integral that we need to evaluate is given by 

^{t)= I Y{s)Y{t-s)Y'^(s)ds. (B.9) 

Jo 

Using the periodicity of the matrix T(t - s), this integral can be 
recast as a finite one over the period of the limit cycle: 



O, 



1 nT„ 

= Yijit - s)e^^<^-^J^'ds 

1 - PiPj Jo 



(B.IO) 



Therefore, the final expression for the correlation matrix is 



c(t + T,t)^ p(t + T)Y(T)<:>(t)P(ty 



(B.ll) 



So we can obtain the correlation matrix as an integral, but this 
has to be evaluated numerically because the neither the limit- 
cycle solutions nor P(t) can be obtained in closed form. 

References 

Alonso, D., McKane, A.J., Pascual, M., 2007. Stochastic amplification in epi- 
demics. J. R. Soc. Interface 4, 575-582. 

Altizer, S., Dobson, A., Hosseini P., Hudson, P., Pascual, M., Rohani, P., 2006. 
Seasonality and the dynamics of infectious diseases. Ecol. Lett. 9, 467^84. 

Anderson, D.F., 2007. A modified next reaction method for simulating chem- 
ical systems with time dependent propensities and delays. J. Chem. Phys. 
127, 214107. 

Anderson, R.M., Grenfell, B.T., May, R.M., 1984. Oscillatory fluctuations in 
the incidence of infectious disease and the impact of vaccination: time seiies 
analysis. J. Hyg. Camb. 93, 587-608. 

Anderson, R.M., May, R.M., 1991. Infectious diseases of humans: dynamics 
and control. Oxford University Press, Oxford. 

Aron, J.L., Schwartz, I.B., 1984. Seasonality and period-doubling bifurcations 
in an epidemic model. J. Theo. Biol. 1 10, 665-679. 

Bartlett, M.S., 1957. Measles periodicity and community size. J. R. Stat. Soc. 
Ser. A 120, 48-70. 

Bartlett, M.S., 1960. Stochastic population models in ecology and epidemiol- 
ogy. Methuen, London. 

Bauch, C.T., Earn, D.J.D., 2003a. Interepidemic intervals in forced and un- 
forced SEIR models. Fields Inst. Comm. 36, 33^4. 

Bauch, C.T., Earn, D.J.D., 2003b. Transients and attractors in epidemics. Proc. 
R. Soc. Lond. B 270, 1573-1578. 

Benton, T.G., 2006. Revealing the ghost in the machine: using spectral analysis 
to understand the influence of noise on population dynamics. Proc. Natl. 
Acad. Sci. USA 103, 18387-18388. 

Black, A.J., McKane, A.J., 2010. Stochasticity in staged models of epidemics: 
quantifying the dynamics of whooping cough. J. R. Soc. Interface 7, 1219- 
1227. 

Boland, R.P, Galla, T., McKane, A.J., 2009. Limit cycles, complex Floquet 
multipliers and intrinsic noise. Phys. Rev. E 79, 051131. 

Bolker, B., Grenfell, B., 1995. Space, persistence and dynamics of measles 
epidemics. Phil. Trans. R. Soc. Lond. B 348, 309-320. 

Bolker, B.M., Grenfell, B.T., 1993. Chaos and biological complexity in measles 
dynamics. Proc. R. Soc. Lond. B 251, 75-81. 

Conlan, A.J.K., Rohani, P, Lloyd, A.L., Keeling, M., Grenfell, B.T., 2009. 
Resolving the impact of waiting time distributions on the persistence of 
measles. J. R. Soc. Interface . 

Coulson, T., Rohani, P., Pascual, M., 2004. Skeletons, noise and population 
growth: the end of an old debate? Trends Ecol. Evol. 19, 359-364. 



10 



Dietz, K., 1976. The incidence of infectious diseases under the influence of 
seasonal fluctuations, in: Lecture notes in Biomathematics. Springer- Verlag. 
volume 11, pp. 1-15. 

Durrett, R., Levin, S.A., 1994. The importance of being discrete and spatial. 
Theo. Pop. Biol. 46, 363-394. 

Earn, D., Rohani, R, Bolker, B., Grenfell, B., 2000. A simple model for com- 
plex dynamical transitions in epidemics. Science 287, 667-670. 

Engbert, R., Drepper, P., 1994. Chance and chaos in population biology - mod- 
els of recurrent epidemics and food chain dynamics. Chaos, Solitons & 
Fractals 4, 1147-1169. 

Fairgrieve, T.F., Jepson, A.D., 1991. O.K. Floquet multipliers. SIAM J. Numer. 
Anal. 28, 1446-1462. 

Ferguson, N.M., Anderson, R.M., Gamett, G.R, 1996. Mass vaccination to 
control chickenpox: The influence of zoster Proc. Natl. Adac. Sci. USA 93, 
7231-7235. 

Forgoston, E., Billings, L., Schwartz, LB., 2009. Accurate noise projection for 
reduced stochastic epidemic models. Chaos 19, 043110. 

Gardiner, C.W., 2003. Handbook of stochastic methods. Springer 3rd edition. 

Gillespie, D., 1976. A general method for numerically simulating the stochastic 
time evolution of coupled chemical reactions. J. Comp. Rhys. 22, 403^34. 

Greenman, J.V., Benton, T.G., 2005. The impact of environmental fluctuations 
on structured discrete time population models: resonance, synchrony and 
threshold behaviour. Theo. Pop. Biol. 68, 217-235. 

Grenfell, B., Bjomstad, O., Finkenstadt, B., 2002. Dynamics of measles epi- 
demics: scaling noise, determinism and predictability with the TSIR model. 
Eco. Monographs 72, 185-202. 

Grenfell, B., Bjomstad, O.N., Keppey, J., 2001. Travelling waves and spatial 
hierarchies in measles epidemics. Nature 414, 716-723. 

Grenfell, B., Harwood, J., 1997. (Meta)population dynamics of infectious dis- 
eases. Trends Ecol. Evol. 12, 395-399. 

Griffiths, D., 1973. The eft'ect of measles vaccination on the incidence of 
measles in the community. J. R. Statist. Soc. A 136, 441. 

Grima, R., 2009. Noise-induced breakdown of the Michaelis-Menten equation 
in steady-state conditions. Rhys. Rev. Lett. 102, 218103. 

Grimshaw, R., 1990. Nonlinear ordinary dift'erential equations. Blackwell, 
Oxford. 1 St edition. 

Hagenaars, T, Donnelly, C, Ferguson, N., 2004. Spatial heterogeneity and the 
persistence of infectious diseases. J. Theo. Biol. 229, 349-359. 

He, D., Earn, D.J.D., 2007. Epidemiological eftects of seasonal oscillations in 
birth rates. Theo. Pop. Biol. 72, 274-291. 

van Kampen, N.G., 1992. Stochastic processes in physics and chemistry. Else- 
vier, Amsterdam. 

Keeling, M., 2000. Metapopulation moments: coupling, stochasticity and per- 
sistence. J. Ani. Ecol. 69, 725-736. 

Keeling, M., Rohani, P., 2002. Estimating spatial coupling in epidemiological 
systems: a mechanistic approach. Ecol. Lett. 5, 20-29. 

Keeling, M., Rohani, P., 2007. Modelling infectious diseases in humans and 
animals. Princeton University Press, Princeton. 

Keeling, M., Rohani, P., Grenfell, B., 2001. Seasonally forced dynamics ex- 
plored as switching between attractors. Physica D 148, 317-335. 

Kravtsov, Y.A., Surovyatkina, E.D., 2003. Nonlinear saturation of prebifurca- 
tion noise amplification. Phys. Lett. A 319, 348-351. 

Kuznetsov, Y.A., 2004. Elements of applied bifurcation theory. Springer, New 
York. 3rd edition. 

Kuznetsov, Y.A., Piccardi, C, 1994. Bifurcation analysis of periodic SEIR and 
SIR epidemic models. J. Math. Biol. 32, 109-121. 

Lloyd, A.L., Sattenspiel, L., 2009. Spatiotemporal dynamics of measles: syn- 
chrony and persistence in a disease metapopulation, in: Spatial Ecology. 
CRCPress, pp. 251-272. 

London, W.P, York, J. A., 1973. Recurrent outbreaks of measles, chickenpox 
and mumps L Seasonal variation in contact rates. Am. J. Epidemiology 98, 
453^68. 

Lust, K., 2001. Improved numerical Floquet multipliers. Int. J. Bifurcation and 
Chaos 11,2389-2410. 

McKane, A.J., Newman, T.J., 2005. Predator-prey cycles from resonant ampli- 
fication of demographic stochasticity. Phys. Rev. Lett. 94, 218102. 

Nasell, I., 1999. On the time to extinction in recurrent epdemics. J. R. Stat. 
Soc. B 61, 309-330. 

Nasell, 1., 2002. Measles outbreaks are not chaotic, in: Mathematical ap- 
proaches for emerging and reemerging infectious diseases: models, methods 
and theory. Springer- Verlag, pp. 85-114. 



Nguyen, H.T.H., Rohani, P., 2008. Noise, nonlinearity and seasonality: the 
epidemics of whooping cough revisited. J. R. Soc. Interface 5, 403^13. 

Nisbet, R., Gurney, W., 1982. Modelling fluctuating populations. Wiley, New 
York. 

Olsen, L.F., Truly, G.L., Schafter, W.M., 1988. Oscillations and chaos in epi- 
demics: a nonlinear dynamic study of six childhood diseases in Copenhagen, 
Denmark. Theo. Pop. Bio. 33, 344-370. 

Priestley, M.B., 1982. Spectral analysis and time series. Academic Press, Lon- 
don. 

Rand, D.A., Wilson, H.B., 1991. Chaotic stochasticity - a ubiquitous source of 
unpredictability in epidemics. Proc. R. Soc. Lond. B 246, 179-184. 

Rohani, P., Earn, D.J.D., Grenfell, B.T., 1999. Opposite patterns of synchrony 
in sympatric disease metapopulations. Science 286, 968-97 1 . 

Rohani, P., Keeling, M.J., Grenfell, B.T, 2002. The interplay between deter- 
minism and stochasticity in childhood diseases. Am. Nat. 159, 469^81. 

Schenzle, D., 1984. An age-structured model of pre- and post-vaccination 
measles transmission. Math. Med. Biol. 1, 169-191. 

Schwartz, I.B., 1985. Multiple stable recurrent outbreaks and predictability in 
seasonally forced nonlinear epidemic models. J. Math. Biol. 21, 347-361. 

Schwartz, LB., Smith, H.L., 1983. Infinite subharmonic bifurcation in an SEIR 
epidemic model. J. Math. Biol. 18, 233-253. 

Simoes, M., Telo da Gama, M.M., Nunes, A., 2008. Stochastic fluctuations in 
epidemics on networks. J. R. Soc. Interface 5, 555-566. 

Wisenfeld, K., 1985a. Noisy precursors of nonlinear instabilities. J. Stat. Phys. 
38, 1071-1097. 

Wisenfeld, K., 1985b. Virtual Hopf phenomenon: a new precursor of period- 
doubling bifurcations. Phys. Rev A 32, 1744-1751. 

Wolfram Research, 2008. Mathematica. Wolfram Research Inc., Champaign, 
Illinois. Version 7. 

Xia, Y, Bjomstad, O.N., Grenfell, B.T, 2004. Measles metapopulation dynam- 
ics: a gravity model for epidemiological coupling and dynamics. Am. Nat. 
164, 267-281. 



11 



